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Abstract. Neutrino interactions beyond the standard model of particle physics may affect 
the cosmological evolution and can be constrained through observations. We consider the 
possibility that neutrinos possess secret scalar or pseudoscalar interactions mediated by the 
Nambu-Goldstone boson of a still unknown spontaneously broken global U{1) symmetry, as 
in, e.g., Majoron models. In such scenarios, neutrinos still decouple at T ~ IMeV, but 
become tightly coupled again (“recouple”) at later stages of the cosmological evolution. We 
use available observations of the cosmic microwave background (GMB) anisotropies, including 
Planck 2013 and the joint BICEP2/Planck 2015 data, to derive constraints on the quantity 
itv: parameterizing the neutrino collision rate due to scalar or pseudoscalar interactions. 
We consider both a minimal extension of the standard ACDM model, and more complicated 
scenarios with extra relativistic degrees of freedom or non-vanishing tensor amplitude. For a 
wide range of dataset and model combinations, we find a typical constraint < 0.9 x 10“^^ 
(95% G.L.), implying an upper limit on the redshift of neutrino recoupling < 8500, 
leaving open the possibility that the latter occured well before hydrogen recombination. In 
the framework of Majoron models, the upper limit on roughly translates on a constraint 
g < 8.2 X 10“^ on the Majoron-neutrino coupling constant g. In general, the data show a 
weak (~ Icr) but intriguing preference for non-zero values of with best hts in the range 
itu — (0.15 — 0.35) X 10“^’^, depending on the particular dataset. This is more evident 
when either high-resolution GMB observations from the ACT and SPT experiments are 
included, or the possibility of non-vanishing tensor modes is considered. In particular, for 
the minimal model KCD^l+^yy and including the Planck 2013, ACT and SPT data, we 
report y^^, = (O.ddlggg) x 10“^^ (300 < Zy^ec ^ 5500) at 68% conhdence level. 
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1 Introduction 

Cosmological observations are a powerful probe of neutrino physics. To date, all the avail¬ 
able cosmological data are consistent with the expectation, based on the standard model 
of particle physics, that the Universe is filled with a thermal background of relic neutrinos 
with Tjy = 1.9 K, belonging to three families, each contributing 113 particles per cm^ to the 
total neutrino abundance. In the early Universe, neutrinos are kept in thermal equilibrium 
with the cosmological fluid by weak interactions until the expansion brings the temperature 
down to T ~ 1 MeV. At later times, the interaction rate becomes too small and neutrinos 
decouple. However, since the decoupling happens when neutrinos are ultrarelativistic, they 
keep a thermal spectrum whose temperature scales as the inverse of the cosmological scale 
factor a. Thus, in the standard cosmological model (SCM), the only free parameters the 
‘neutrino sector’ of the model are the masses of the three eigenstates. Oscillation exper¬ 
iments have measured mass differences, showing that at least two of the eigenstates have 
non-vanishing masses. However, both the absolute scale and the hierarchy of the masses 
remain unknown. Present-day cosmological observations can tightly constrain the sum of 
the masses [1, 2]: Planck measurements of the temperature and polarization anisotropies of 
the cosmic microwave background (CMB), when combined with other astrophysical datasets, 
imply [3]. 

This simple picture, other than being theoretically well-grounded, is perfectly consistent 
with all the available data. It is however desirable to keep an open mind and test more 
complicated scenarios for the neutrino sector of the SCM. For example, cosmological data 
can be used, among others, to constrain lepton asymmetry [4], possible deviations of the 
neutrino spectrum from equilibrium [5], and also the properties of a sterile neutrino with ~ 
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eV mass [6-20], that could possibly explain reactor anomalies (see e.g. Refs. [21, 22] for a 
review). In this paper, we will consider the possibility that neutrinos have interactions beyond 
the standard model of particle physics (that for simplicity we shall call “hidden” or “secret” 
interactions) and study the constraining power of cosmological observations with respect 
to such a scenario. In particular, we will consider a specific version of secret interactions, 
that is however representative of a large class of models, namely a (pseudo)scalar interaction 
mediated by the Nambu-Goldstone boson of a hitherto unknown broken U{1) symmetry, like 
in Majoron models [23-25]. These models are very well-motivated from the point of view 
of particle physics and have the peculiar feature that the ratio between the interaction and 
expansion rates increases with time, so that neutrinos may actually become kinematically 
coupled again in the late stages of the cosmological evolution. 

Non-standard neutrino interactions have first been considered in a cosmological context 
in Ref. [26], discussing the possibility that neutrino decay induced by the new interaction 
would lead to a neutrinoless Universe. Limits on neutrino-neutrino scattering induced by 
non-standard interactions (either Majoron-like, as those considered in this paper, or Fermi- 
like) from cosmological observations have been derived in Refs. [27-29] and, more recently, 
in Ref. [30], using data from the Planck 2013 release. In these papers, the neutrino fluid 
is modeled as abruptly changing from collisionless to perfectly tightly coupled (or viceversa 
in the case of Fermi-like interactions) at a given transition redshift, that represents the 
parameter actually constrained by the data. A complementary approach consists in deriving 
limits on phenomenological quantities parameterizing the effective sound speed and viscosity 
of the neutrino fluid [31-37]. As noted by a few authors, however, this approach does not 
always provide an accurate representation of the collisional regime [29, 59]. For this reason 
we avoid resorting to it throughout this paper. Instead, we derive limits on the strength of 
neutrino non-standard interactions by directly modifying the Boltzmann equation in order 
to account for neutrino collisions, without assuming a sudden transition between the two 
limiting regimes (free-streaming and tight coupling). We also consider extended models 
allowing for extra relativistic species or tensor pertubations. 

This paper is structured as follows. In Sec. 2 we briefly introduce the theoretical frame¬ 
work that describes the hidden interactions of interest. In Sec. 3 we review the Boltzmann 
formalism for interacting neutrinos. In Sec. 4 we describe the method used to compute 
the impact of interacting neutrinos on the evolution of cosmological perturbations and on 
the CMB observables, and to derive constraints on the strength of the interaction, that are 
discussed in Sec. 5. Finally, in Sec 6 we draw our conclusions. 

2 Hidden neutrino interactions 

We consider neutrinos interacting with a light boson (p through simple scalar hij and pseu¬ 
doscalar gij couplings, as described by the following Lagrangian [23-25]: 

£ = hijViVjp + gijva^vjp + h.c. , (2.1) 

where the indices i, j run over the neutrino mass eigenstates. This kind of interaction 
allows for the binary processes shown in Fig. I, i.e. v + v p + cp (neutrino annihilation 
to (/)’s), u + p V + p (neutrino-()) scattering), u + v u + u (neutrino-neutrino scattering 
mediated by a scalar boson exchange), as well as for neutrino decay u ^ v' + p. 

Neutrino scalar and pseudoscalar couplings are constrained by laboratory searches for 
neutrinoless double beta decay (Oz^/3/3), and by supernovae observations. For example, in 
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Figure 1. Feynman diagrams for the binary processes allowed by the Lagrangian (2.1). Time goes 
from left to right. From left to right: v-v scattering (s and t channels), v-(\) scattering, vi> annihilation 
to (/)’S. 


addition to the simplest Oi//3/3 decay mode 

(A, Z) ^ (A ^ + 2) + 2e-, (2.2) 

whose existence only requires the neutrino to be a Majorana particle [38], modes in which 
one or two additional (p bosons are emitted 

(^, Z) —> [A, Z + 2) + 2e (p , (2-3) 

[A, Z) — > [A, Z + 2) + 2e + 2</>, (2-4) 

are possible if neutrinos possess (pseudo)scalar couplings. Oz//3/3 experiments yield constraints 
on the effective ^-neutrino coupling constant {gee) < (0.8 — 1.6) x 10“®, depending on the 
theoretical model [39, 40]. The quantity gee is the e — e entry of the coupling matrix in the 
weak base, related to the couplings gij in the mass basis through the elements of the neutrino 
mixing matrix. 

Neutrino decays v ^ v' + p can also be important in the high-density supernova envi¬ 
ronment [41-44]. In the case of Majoron models, limits on Majoron-neutrino couplings from 
observations of SN 1987A were derived in Ref. [42]. It has been shown there that p emission 
would shorten too much the observed neutrino signal from SN 1987A if 3 x 10“^ ^ 5 ^ 2x 10“® 
(here g denotes the largest element of the coupling matrix g^iB in the weak base), thereby 
excluding this region. Moreover, the observed Ue flux from SN 1987A can also be used to fur¬ 
ther constraint gn < 10“^. These limits, together with those from Oz^/3/3 decay experiments 
available at that time, were combined and translated into the mass basis in Ref. [43]. 

Scalar and pseudoscalar neutrino couplings can also be relevant in a cosmological con¬ 
text, since collisional processes induced by the new interaction would affect the evolution of 
perturbations in the cosmological neutrino fluid. In general, the cross section for a binary 
process mediated by a massless boson has the form fjbin ~ g'^/s in the ultrarelativistic limit 
(apart from numerical factors) with g being the largest value of the Yukawa matrix (we do 
not distinguish between scalar and pseudoscalar couplings in the following), and y/s is the 
center-of-mass energy. Thus, in thermal equilibrium, the rate for a binary process is 

Tbin — (u'binU)rieq OC. g T , (2-5) 

since the equilibrium neutrino abundance Ueq oc T^, and s ^ . 

Interactions are of cosmological significance when the ratio T /H between the interaction 
and Hubble expansion rates is of order unity or larger. The expansion rate scales as H r2 
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(H J'3/2) the radiation- (matter-) dominated era. The fact that the interaction rate 

(2.5) scales like T has interesting phenomenological consequences: since the ratio Tbin/^^ 
actually increases with decreasing temperature, neutrinos, having decoupled in the early 
Universe, could recouple at late times due to scattering/annihilation processes mediated by 
4 >. 

Let us describe in more detail what happens once neutrino secret interactions be¬ 
come effective. We define the recoupling redshift of neutrinos ^^recj through the condition 
rbin(^!^rec) = H{zjjrec)- At z < Zjyrec, tbe neutrino free-streaming length rapidly drops well 
below the Hubble length due to scatterings. Hence, the neutrino contribution to the cosmic 
shear becomes negligible. This effect should be observable in the CMB anisotropy spectrum, 
provided recoupling happens close enough to recombination. Once recoupling becomes ef¬ 
fective, one should expect production of 0’s through vD annihilation, and their subsequent 
thermalization through scatterings^. This allows us to describe the tightly coupled z^’s 
and 0’s as a single fluid. Note that if late 0 production, as we assume, happens when both 
neutrinos and 0’s are in the ultrarelativistic regime, the total energy stored in relativistic 
species does not change, so A^eff remains constant. This holds as long as chemical equilibrium 
is mantained. Once the temperature of the fluid falls below the neutrino mass, annihilation 
processes start to deplete the neutrino abundance, and the energy stored in the neutrino rest 
mass would end up increasing Neff [45, 46]. We neglect this effect in what follows, assuming 
that neutrino masses are small enough to become cosmologically significant well after recou¬ 
pling. Finally, we also neglect the possibility of neutrino decay, implicitly assuming that the 
off-diagonal couplings are small. 

3 The Bolzmann equation 
3.1 Formalism 

The phase-space evolution of the components of the primordial plasma is described by the 
Boltzmann equation, that can be generically written as: 

(3.1) 

where / is the phase-space distribution function (DF), and C[f] is the collision operator, 
describing interactions between particles, and r is a generic time variable, that we here 
take to be conformal time. In fact, the cosmological evolution is described by a set of 
Boltzmann equations, one for each component of the cosmological fluid, coupled between 
them by gravity [hidden in the LHS of Eq. (3.1)], and possibly by the collision term (as in 
the case of baryons and photons, coupled by Thomson scattering). This system of equations 
must be complemented by the Einstein equations describing the evolution of metric variables. 

In a perturbed Eriedmann-Robertson-Walker Universe, the DE can be conveniently 
written as the sum of an unperturbed part /o (that does not depend on position nor on the 
direction of momentum, due to the homogeneity and isotropy of the background spacetime) 
plus a small perturbation 6f = /odz: 

f{x,q,T) = fo{q)[l + 4'(f,g,T)] , (3.2) 

^We neglect the population of 0 bosons produced during reheating and the correspoding contribution 
to Neff, since their energy density would be diluted by the entropy produced by standard model particles 
annihilations. 


Dl 

Dt 


= C[f] 



where we are using comoving coordinates x* and momenta qi = qrii {rii being a unit vector) 
as the spatial and momentum coordinates respectively. 

In the synchronous gauge, the perturbed metric is written in the form 

ds‘^ = a^(r) [—+ {dij + hij) dx^dx^] , (3.3) 

where a(r) is the cosmic scale factor and hij are the metric perturbations. With this choice 
of the gauge, the perturbed Boltzmann equation in /c-space takes the form (for the sake of 
simplicity, we use the same symbol for 'k and its Fourier transform): 


oik q 

— -h ikfi-^ + 

OT e 


ding 


h + Qfj 2 




(3.4) 


where e = \/q^~-\-~a?m^ {m being the mass of the particle), = k ■ h is the angle between 
the perturbation wave number and the particle momentum, and h and r] are the scalar 
components of the metric perturbation hij in /c-space: 


hij{x, t) 



kikjh{k, t) + {kikj 


1 

3 


dij) 6r]{k, t) 


(3.5) 


For the moment, let us consider a collisionless fluid, i.e., C[f] = 0. In the case of massless 
particles, e = q and the Boltzmann equation (3.4) can be further simplified by integrating 
out the momentum dependence of the DF. Defining 


F{k, h, t) = 


/gVo(g)4'(fc, q, t) dq 

Jq^foiq) dq 


(3.6) 


and further expanding the angular dependence of F" in a series of Legendre polynomials: 


F{k, h, t) = + l)F,{k, T)Pe{fi) (3.7) 

e=o 


the following Boltzmann hierarchy is obtained: 


4 2 ■ 

d = — e--h. 
3 3 

0 = k^ l-S-U 


4 3 2-4 

n = —6 - kF^ H- h H—h , 

15 10 '^15 5'’ 


Fp = 


k 


2£+l 


iFp_^ -ii + l)F,+i ii>3), 


(3.8a) 

(3.8b) 

(3.8c) 

(3.8d) 


where d = Fq, 6 = (3/4)/cFi and 11 = F 2/2 (our 11 corresponds to a in the notation of Ref. 
[49]). 


3.2 Boltzmann equation for interacting neutrinos 

In order to study the behaviour of non-standard interacting neutrinos, we need to specify 
the form of the collision term on the RHS side of the Boltzmann equation. In principle, 
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we should compute the collision integrals given a specific form of the interaction lagrangian. 
Detailed calcnlations of the collision term for neutrino-neutrino interactions mediated by a 
scalar particle have been presented, under some assumptions, in Ref. [59]. In that work, 
the collision integrals are simplified by redncing their dimensionality, but yet they cannot 
be computed analytically, leaving the Boltzmann equation in its integro-differential form. 
Moreover, the presence of the collision terms prevents from integrating out the momentum 
dependence in the Boltzmann equation also in the massless case, increasing the numerical 
complexity of the problem. 

In this work, instead, we will pursue a simpler approach and use the relaxation time 
approximation as done in [60]. This approximation amounts in writing the collision integral 
as a damping term, proportional to the inverse mean free time between collisions Tc, i.e. 

C[/](x--<5/. (3.9) 

'^C 

where = aT = an{av). Here T = n{av) is the interaction rate in the comoving frame, and 
the scale factor comes from the fact that we choose conformal time as the time coordinate. 

Following the same steps that led from Eq. (3.4) to Eqs. (3.8), it is straighforward to 
show that a collision term of the form (3.9) would result in additional terms proportional 
to —aTFi on the RHS of the equation for Fi. Notice, however, that conservation of the 
nnmber of particles (as in a 2 -f-)- 2 scattering process) and conservation of momentum imply 
respectively 


j q^C[f] dndq = 0, 
f nq^C[f] dQ dq = 0 , 


(3.10) 

(3.11) 


i.e., no collision terms should appear in the monopole and dipole evolution equations. Eor 
this reason we set to zero the collision terms for i = 0, 1 , so that the Boltzman hierarchy for 
interacting neutrinos reads [we choose the proportionality constant in Eq. (3.9) to be unity] 


4 2 ■ 

,5 = — e--h, 

3 3 

e = k^ (^^6-u 

4 3 2-4 

n = —0 - kF^ H h H— f] — aFH , 

15 10 15 5 ' 


Fe = 


k 


2£ + l 


^F,_i-(£+l)F>+i -aTFi (£ > 3) 


(3.12a) 

(3.12b) 

(3.12c) 

(3.12d) 


We have noted above that, once secret interactions become cosmologically significant, 
neutrino-neutrino annihilations start to prodnce (p bosons that are rapidly brought into ther¬ 
mal equilibrium with neutrinos themselves by scatterings. This allows to treat the neutrino-(/> 
system as a single, tightly coupled self-interacting fluid, whose perturbations still evolve ac¬ 
cording to Eqs. (3.12). The conservation of particle number and momentum still holds, at 
the level of the coupled fluid, justifying the vanishing I" = 0, 1 collision terms also in this 
regime. 
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4 Method 


In the following, we parameterize the strength of the non-standard interaction by generically 
writing the coefficient of the collisional damping terms appearing in the RHS of Eqs. (3.12) 
for £ > 2 as 

= (4.1) 

This has the right energy dependence for an interaction mediated by a light (pseudo)scalar 
boson in the ultrarelativistic limit. With such a choice, the parameter is roughly pro¬ 
portional, through a dimensionless factor, to the largest Yukawa coupling appearing in the 
Lagrangian (2.1). The exact relation between and the gij's depends on the details of the 
underlying particle physics model. Nevertheless, in the following, we shall loosely refer to 
7 j/i/ as the coupling constant for non-standard neutrino interactions. 

We use a modified version of the publicly available Boltzmann code camb [48] in order 
to solve the perturbation evolution and compute the CMB power spectra for given values 
of the cosmological parameters. In particular, we have modified the evolution equations for 
interacting massless neutrinos described in Sec. 3.2 in order to include a collision term with 
the form (4.1). 

4.1 Cosmological data analysis 

The baseline model considered in this paper is described by the following parameter set: 

{Ulb, Uc, 9, Tree, Us, log[10^°4s], 7 ]^^} , (4.2) 

where ojb and lOc are the physical baryon and cold dark matter density, respectively, 9 is the 
angle subtended by the sound horizon at the time of recombination. Tree is the optical depth 
to reionization, Ug and Ag are the spectral index and amplitude of the primordial spectrum 
of scalar fluctuations (both evaluated at the pivot wavenumber ko = 0.05 Mpc“^), and 7 ^^, 
parameterizes the strength of non-standard neutrino interactions, as per Eq. (4.1). These are 
the parameters that are varied in the MC run and that implicitly take flat priors. We assume 
flat geometry, massless neutrinos and adiabatic initial conditions. In the baseline model, we 
also fix Yeff = 3.046 and ignore tensor modes by setting the tensor-to-scalar ratio r = 0. We 
refer to this model as h.CDM.+^,yy. When is fixed to zero, it reduces to the standard 
ACDM model with weakly-interacting neutrinos. We also consider extensions to the baseline 
model, by allowing Yeg or r to vary, referred to as ACDM-|- 7 ,yjy-|-Yeg and ACDM-|- 7 jyi,-|-r, 
respectively. In the following we shall also presents results for derived parameters of in¬ 
terest, like the present Hubble parameter Yq, the tensor-to-scalar ratio ro .002 evaluated at 
0.002 Mpc“^ and the neutrino recoupling redshift Zy^ec- la Tab. 1 we summarize our choice 
of parameterization. In addition to the parameters listed in the table, we also vary a number 
of “nuisance” parameters describing residual foregrounds and instrumental characteristics. 

In order to derive constraints for the parameters of the model, we perform a Markov 
Chain Monte Carlo (MCMC) analysis through the publicly available CosmoMC code [50], inter¬ 
faced with our modified version of camb. CosmoMC uses a modihed version of the Metropolis- 
Hustings algorithm [51] to sample the full posterior distribution of the parameters given 
the data (described in the following section). Lower-dimensional posterior distributions are 
obtained by marginalizing over unwanted parameters. When quoting credible intervals, we 
apply the following rule: if both edges of the 95% credible interval are distinct from the prior 
edges, we quote constraints in the form mean ± 68 % uncertainty; on the contrary we quote 
the 95% upper or lower limit. 
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Table 1. Cosmological parameter used in the analysis. The upper part of the table lists the base 
parameters, i.e., those with uniform priors that are varied in the Monte Carlo runs. In the baseline 
model we only consider the first seven of those listed here; A^eff and r are varied in extensions of 
this model. The lower part lists derived parameters of interest, for which we also compute credible 
intervals. For each parameter, we quote the initial prior range (for base parameters only). 


Parameter 

Definition 

Prior range 


Present baryon density 

[0.005, 0.1] 


Present dark matter density 

[0.005, 0.1] 

100 e 

Angular size of the sound horizon at recombination 

[0.5, 10] 

Tree 

Optical depth to recombination 

[0.01, 0.8] 

Us 

Spectral index of scalar perturbations 

[0.5, 1.5] 

ln(10^%) 

Log amplitude of scalar perturbations at ko = 0.05 Mpc“^ 

[2.7, 4.0] 


Strength of non-standard neutrino interactions^ 

[0, 2] 


Effective number of neutrino families^ 

[0, 5] 

r 

Tensor-to-scalar ratio at ko = 0.05 Mpc“^ ^ 

[0, 2] 

Ho 

Hubble constant at a{t) = 1 

_ 

ro.oo 2 

Tensor-to-scalar ratio at ko = 0.002 Mpc“^ 

- 

^urec 

Redshift of neutrino recoupling 

- 


4.2 Cosmological data 

We consider the data on CMB temperature anisotropies released by the Planck satellite in 
2013 [52, 54], supplemented by the 9-year polarization data from WMAP [47], as well as 
additional temperature data from high-resolution CMB experiments, namely the Atacama 
Cosmology Telescope (ACT) [55] and the South Pole Telescope (SPT) [56]. The purpose 
of considering the ACT and SPT data is mainly to improve constraints on the unresolved 
foregrounds. 

The likelihood functions associated to these datasets are evaluated and combined using 
the likelihood code distributed by the Planck collaboration, described in Ref. [54], and 
publicly available at Planck Legacy Archive^. This likelihood uses Planck TT data up to a 
maximum multipole number of ^max = 2500, and WMAP 9-year polarization data (WP) up 
to £ = 23 [47], as well as ACT data in the range 1000 < £ < 9440 [55] and SPT data in the 
range 2000 <i< 10000 [56]. 

We also use the likelihood recently released by the joint Planck and BICEP2/Keck 
effort [61]. This likelihood is based on CMB polarization observations from the BICEP2 
field and uses corresponding data from Planck 2015 217 and 353 GHz channels to account for 
contamination from polarized Galactic dust. It is limited to the multipole range 20 < i < 200. 

®http://pla.esac.esa.int/pla/ 
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5 Results 


5.1 Impact on the evolution of pertubations 

Once neutrinos become interacting again, their anisotropic stress (as well as all the higher- 
order moments of the distribution function) is suppressed and eventually vanishes once the 
tightly coupled regime {Tuv/H —)> oo) is reached. At the same time, the damping of density 
perturbation caused by neutrino free streaming is no more effective and the neutrino fluid 
undergoes acoustic oscillations. In order to illustrate this, in Fig. 2 we show the evolution of 
density (Jj/) and shear (11,^) perturbations for three different wave numbers {k = 5x 10“^, 5 x 
10- 5 X 10 ^ Mpc ^), in the case of non-interacting neutrinos {i.e., 'ji.i, = 0), and for two 

finite values of the secret interaction strength, namely = 1.2 x 10“^, 2 x 10“'^. The latter 
values correspond to a redshift of neutrino decoupling z^^ec — 1300 and 1.8 x 10^, respectively 
(in this subsection we take the following benchmark for the cosmological parameters: ojb = 
0.0227, LOc = 0.124, h = 0.68, Tree = 0.09, n* = 0.96, Ag = 2.1 x 10^®, Agff = 3.046, 

r = 0). The upper panels of Fig. 2 show perturbation evolution for the largest scale, 

k = 5 X 10“^Mpc“^, entering the horizon around the time of hydrogen recombination at 
z ~ 1100. For the larger value of the coupling constant, when the mode enters the horizon 
neutrinos are already completely recoupled, and shear oscillations are overdamped. For 
= 2 X 10“^, on the other hand, neutrinos are only partially coupled at the time of horizon 
crossing and this results in a significant but not complete suppression of shear perturbations 
with respect to the non-interacting case. The evolution of density perturbations mirrors that 
of the shear: when the dissipation normally associated to neutrino free-streaming is absent, 
undamped acoustic oscillations set on in the fluid, so that density perturbations are actually 
boosted by increasing In the middle row, we show results for k = 5 x 10“^Mpc“^. 
This scale crosses the horizon at z ~ 2 x 10^. In this case, the shear is again suppressed 
as soon as the mode enters the horizon for = 2 x 10“^, since neutrinos are recoupling 

roughly at the same time, while for = 1.2 x 10“^ the effect of interactions only kicks 

off after perturbations have been damped by the expansion. Finally, for the smallest scale 
under consideration {k = 0.5Mpc“^), shown in the bottom row and entering the horizon at 
z ~ 2 X 10^, neutrino recoupling happens in both cases when the perturbation is well inside 
the horizon, making the evolution very similar to the non-interacting case. 

In Fig. 3, we show how secret neutrino interactions affect the CMB power spectra. 
We have computed the power spectra for the same models shown in Fig. 2, i.e. ji/i, = 
0, 1.2x10’ -7, 2x 10 The effect of neutrino interactions on the temperature power spectrum 
shown in the upper panel of Fig. 3 is to boost the overall amplitude of the spectrum as 
is increased. The increased power is caused by the increased density fluctuations due to 
the absence of neutrino free streaming. A similar behaviour is observed in the TE power 
spectrum, as shown in the lower panel of the same figure. 

5.2 Constraints on cosmological parameters 

We are now ready to present the constraints that CMB data provide on the non-standard 
coupling constant juu. In the following, we will quote 68% CL uncertainties, unless we are 
dealing with upper limits, in which case we quote 95% credible intervals. The results shown 
in the following are summarized in Tabs. 2 and 3. 

Let us start by considering the simplest extension of the standard cosmological model, 
labelled ACDM -|- Considering only Planck temperature and WMAP polarization data 
(Planck-|-WP), we obtain < 0.85 x 10“^^, or equivalently juu < 1-71 x 10“’^, which implies 
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Figure 2. Evolution of neutrino density (left column) and shear (right column) perturbations for 
three different modes k = 0.005, 0.05, 0.5Mpc“^. We compare the evolution for standard, weakly- 
interacting neutrinos (blue solid curves) with the case of non-standard neutrino interactions with 
= 1.2 X 10“^ and = 2 x 10“^ (green dashed and red dotted curves respectively). Non¬ 
standard interactions suppress shear perturbations while boosting density perturbations after the 
time of neutrino recoupling (see text for discussion). 


a neutrino-neutrino recoupling at Zuiec ^ 8800 (here and for the rest of the section, we will 
fix the other parameters to their best estimates when translating limits on to limits on 
Zi/rec)- Adding the ACT and SPT datasets (“highL”) shifts the distribution to larger values 
of the coupling constant, yielding 7 ^^, < 0.96 x 10“^’^ {'y^u < 1-76 x 10“^ and Zuiec < 10^). 
In the upper panel of Fig. 4 we show the posterior distributions for 'y^^ in the ACDM-l-yj^i/ 
model, for the Planck-bWP and Planck-|-WP-|-highL datasets. Both posteriors are quite 
asymmetric and have a peak at non-zero values of the coupling constant, respectively at 
7 ^^, = 0.24 X 10“^^ and 0.36 x 10“^^, corresponding to z^rec — 2800 and 1700. Interestingly 
enough, there is a weak (at the ~ lex level) preference for non-zero values of 7 ^^,: at 68 % 
CL, we find 7 ^^ = (O.SdD^g;^^?) x lO-^^ (Planck+WP) and 7 ^^ = (0.4441°;^®^) x lO'^^ 
(Planck-|-WP-|-HL). The 68 % lower limit in the latter case corresponds to Zj^rec ^ 300. 

We have also constrained the number of relativistic species in conjuction with 7 ,^ 1 ,. In the 
framework of this ACDM-|- 7 jyi/-|-Aefr model, we find a 95% credible interval < 0.75 x 10“^'^, 
or 'y^u < 1.65 x 10“^ from Planck-|-WP. This value provides a neutrino-neutrino recoupling 
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Figure 3. CMB temperature (left panel) and temperature-£1 polarization (right panel) power spectra 
for neutrinos with non-standard interactions. As in Fig. 2 we show results for = {0, 1.2 x 10“^, 2 x 
10 -^}. 


at Zjyrec ^ 7400. Also in this case, the addition of the ACT and SPT datasets weakens the 
constraints on the coupling constant, yielding 7 ^^, < 0.93 x 10“^”^ {'jiyi, < 1.75 x 10“’^ and 
Zuiec ^ 9800) at 95% CL. For what concerns the effective number of relativistic species, we 
find A^eff = 3.44l°;37 (Planck+WP) and A^eff = 3.27l[J;| (Planck+WP+HL). This is very 
much consistent with the corresponding values found by the Planck collaboration in the Aeff 
extension of the ACDM model for the same datasets [53]. The posterior distributions for 'y^i, 
in the ACDM-|- 7 i/i/-|-Al^efr model, for the datasets under consideration are presented in the right 
pannel of Fig. 4. The maximum probability is obtained for 7 ^^ = 0.18 x 10“^^ {zurec — 1000) 
and 7 yj, = 0.27 x 10“^^ (z^rec — 2000) for Planck-|-WP and Planck-|-WP-|-HL, respectively. 
The presence of additional relativistic degrees of freedom reduces the preference for non-zero 
values of the coupling constant: in the ACDM-t-yj^j/ -|- A^eff model, 7 ^^, is consistent with 
zero at below the ~ Icr level for Planck-I-WP, while the 68 % interval for Planck-|-WP-|-HL 
is = 0.403lg'3g3, thus shifted to lower values with respect to the corresponding interval 
in the ACDM-l- 71 /j^ model. In Fig. 5, we show the most significant correlations between 7 ^^^ 
and other parameters, namely 9, and A"eff- The correlations with the angle 

9 subtended by the sound horizon at recombination and with the amplitude lO^Age”^"^ are 
particularly evident. We argue that the pattern leading to these correlations is the following: 
the overall amplitude of the spectrum increases for larger values of 7 ^j^, while the position 
of peaks and dips remains unchanged. This can be directly compensated by a lower value 
of Alternatively, increasing (or decreasing A^gff if the model allows), lowers 

the height of the first few peaks but shifts their position to lower multipoles; increasing 9 
moves the peaks back to their original position. To support our reasoning, we show in the 
left panel of Fig. 6 a scatter plot of samples from Planck-|-WP chains in the — 9 plane, 
color-coded by 7 ^^,, compared with the 68 % and 95% confidence region for the ACDM model. 
It is clear from this plot that models with 7 ^^, > 0 can be made in agreement with the data by 
increasing both and 9. The right panel of the same figure shows the same information 
in terms of Hq instead than 9. 

Finally, we have also considered the possibility of a non-vanishing amplitude of tensor 
modes. We label this model as ACDM-|- 7 jyi/-|-r. In this case, in addition to Planck-|-WP, we 
have also used the joint BICEP2/Planck 2015 (BKP) dataset. At 95% CL, we find 7 ^^ < 
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ACDM-h7^^ 

ACDM±7,, ± A^eff 


Planck-|-WP 

Planck±WP 

Planck±WP 

Planck±WP 



±highL 


±highL 

Parameter 






0.02214 ± 0.00029 

0.02220 ± 0.00029 

0.02244 ± 0.00041 

0.02237 ± 0.00040 


0.1218 ± 0.0029 

0.1220 ±0.0029 

0.1263 ±0.0054 

0.1247 ±0.0048 

100 0 

1.04195trooo85 

1.04210t°;°°°™ 

i.U^l^O_0.00097 

1.0418ltroo?o3 

^rec 

0.09lini4 

0.093tHM 

u.uyo_o.oi6 

0.095ini^ 

Us 

0.9641 ± 0.0074 

0.9629 ± 0.0072 

0.979 ±0.016 

0.971 ±0.015 

log[10iOA«] 

O.U<y_Q Q26 

3.080t°;°p 

0 1 (A9+O.O33 
O.lUZ_0 037 

O.uyz_o 036 


< 0.85 

< 0.96 

< 0.75 

< 0.93 

Ks 

3.046 

3.046 

4 44+0-3'^ 

q 9'7+0.33 

Ho [km/sec/Mpc] 

67.4 ± 1.2 

67.4 ± 1.2 



10'^7i,i. 

< 1.71 

< 1.76 

< 1.65 

< 1.75 


Table 2. Constraints on cosmological parameters for the ACDM+7i^,y and ACDM+7^;^+A"eff models 
from the analysis of the Planck+WP and Planck+WP+highL datasets. We quote 68% C.L., except 
for upper bounds, which are 95% C.L. 



Figure 4. One-dimensional posterior distribution for 7^^ obtained from the Planck-|-WP (blue solid) 
and Planck-l-WP-l-highL (red dashed) datasets. The shaded regions denote the 68% credible interval. 
Left panel: ACDM -|- 71,^. Right panel: ACDM-|-7i,i,-|-A^eff ■ 
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Figure 5. 68 % and 95% confidence regions for selected parameter pairs involving 7 ^^ in the 

ACDM+ 7 ,yi, (empty contours) and ACDM+ 7 y;, + A^eff (filled contours), for Planck+WP (blue) and 
Planck+WP+HL (red). 


0.90 X 10“^^ (7j/i/ < 1.73 X 10“^ and ^i/rec ^ 9300) for Planck+WP and 7)^^ < 0.87 x 10“^^ 
{ivv < 1-72 X 10“^ and z^rec ^ 9000) for Planck+WP+BKP. Even in this extension of the 
ACDM model, cosmological data slightly prefer a non-zero value of the coupling: at 68% 
CL, we have 7^^, = (0.408lg 357) x 10“^’^ for Planck+WP and = ( 0 . 395 ^q 342 ) x 10“^^ 
for Planck+WP+BKP. The posterior probability for in the two cases is shown in Fig. 7. 
For what concerns the tensor-to-scalar ratio, we find r < 0.14 and r < 0.10 for Planck+WP 
and Planck+WP+BKP, respectively. Both values are consistent with those reported in Refs. 
[53, 61] for the ACDM+r model. 
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Figure 6. Left panel: Samples from Planck+WP chains in the plane, color-coded by 7 ^,^ for 

the ACDM-|- 7 ;^jy model. The 68 % and 95% credible regions for the same dataset in the ACDM model 
are also shown (solid lines). Notice how larger values of the coupling constant for secret interactions 
require larger values of both and 0, as explained in the text. Right panel: The same as the left 
panel, but in the Q,ch? — Hq plane. 



Figure 7. One-dimensional posterior distribution for obtained from the Planck-|-WP (blue solid) 
and Planck -I- WP-I-BKP (red dashed) datasets, for the ACDM -|- 7^1. -|- r model. The shaded regions 
denote the 68 % credible interval. 


6 Conclusions 

We have derived constraints on non-standard neutrino interactions mediated by a (pseudo)scalar 
massless boson using observations of CMB temperature and polarization anisotropies from 
Planck, WMAP, ACT, SPT and BICEP2/KECK. We have found that, both in a minimal 
extension of the ACDM model and in more complicated scenarios allowing for the presence 
of extra relativistic degrees of freedom or of primordial tensor perturbations, the strength 
of non-standard interactions (expressed through the coefficient of the collision term in the 
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Parameter 


Planck-tWP 

Planck±WP 

±BKP 

100 B 

^rec 

ris 

log[10^°As] 

r 

0.02220 ± 0.00029 

0.1213 ± 0.0029 

1.042111^000^5 

o-09iini4 

0.9669 ± 0.0078 

o n'Tc;+0.025 
O.U<0_Q Q28 

< 0.90 

< 0.14 

0.02217 ± 0.00029 

0.1217 ±0.0029 

1.042031°;^ 

0.0911°:°(3 

0.9658 ± 0.076 

3.078 ±0.025 

< 0.87 

< 0.10 

Hq [km/sec/Mpc] 

67.7 ± 1.2 

67.5 ± 1.2 


< 1.73 

< 1.72 

ro.oo2 

< 0.13 

< 0.09 


Table 3. Constraints on cosmological parameters for the ACDM+ 7 yj^+r model from the analysis 
of the Planck+WP and Planck+WP+BKP datasets. We quote 68% C.L., except for upper bounds, 
which are 95% C.L. 


Boltzmann equation for neutrinos) is constrained at 95% C.L. < 1.7 x 10“^, quite stable 
with respect to the models and datasets considered. This corresponds to a largest redshift 
of neutrino recoupling of Zi/rec — 8500, larger than the value Zy-cec < 1887 found in Ref. 
[30], and shows that the possibility of neutrino recoupling happening before recombination 
is allowed by the data. On the other hand, we confirm the preference, also reported in Ref. 
[30], for non-zero values of the coupling constant. We find best-fit values of in the range 
(1.2 -b 1.4) X 10“^, corresponding to in the range 1300 -b 3300. For comparison, in Ref. 
[30] it is found that the probability distribution peaks in Zy^^c — 1500. In most cases, we find 
"^yy 7 ^ 0 at 68% CL; this tendency is more pronounced when small-scale CMB observations, 
which are sensitive to details of the photon damping regime, are considered, but is alleviated 
in presence of extra relativistic degrees of freedom if one allows for them. On the other 
hand, considering a non-vanishing amplitude of tensor modes, still leads to a preference for 
non-zero coupling at the same level, even for the base Planck-|-WP dataset. 

The exact relationship between our parameter ^yy and the elements of the Yukawa 
matrix gij depends on the details of the underlying particle physics model. As an example, let 
us consider the class of models in which neutrino acquire mass through violation of ungauged 
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lepton number. In this case the neutrino mass eigenstates couple diagonally, to lowest order 
approximation, to the Nambu-Goldstone boson of the broken global symmetry, the Majoron. 
Neutrino masses are proportional to the diagonal couplings: rrii oc ga. Neglecting the small 
off-diagonal couplings, gij = Sijgi, and further assuming that the diagonal ones are of the 
same order of magnitude, gi ~ g, we have that 

Tbin = neq(fTbinw) ~ (1.8 X 10“^) g^T„, (6.1) 

where we have used a = g'^/{32Trs) [62] for the neutrino-neutrino scattering cross section and 
rieq is the abundance of single neutrino family. Comparing this with Eq. (4.1) immediately 
yields g < 8.2 x 10“^. This region partially overlaps with the interval 3 x 10“^ ^ 5 ^ 2 x 10“^ 
excluded by observations of SN1987A [42], although as discussed in Sec. 2 SN observations 
do not directly probe the diagonal elements of the coupling matrix in the mass base. The 
best-fit values for 7 ^^^ translate to <7 ~ (5 -7 6 ) x 10“'^ in the Majoron model, which is also in 
tension with SN bounds, albeit the same remark as above applies. 
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